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ABSTRACT 

In  this  paper,  we  propose  and  validate  a 
computationally  efficient  non-iterative  domain 
decomposition  procedure  for  calculating  bivariate  cubic 
Lx  smoothing  splines.  This  domain  decomposition 
procedure  involves  calculating  local  Lx  smoothing  splines 
individually  on  overlapping  “extended  subdomains”  that 
cover  the  global  domain  and  then  creating  the  global  Lx 
smoothing  spline  by  patching  together  the  local  Lx 
smoothing  splines.  Using  this  procedure,  we  calculate  the 
global  Lx  smoothing  splines  of  one  urban  terrain  data  set 
(Baltimore)  and  one  natural  terrain  data  set  (Killeen, 
Texas).  The  local  Lx  smoothing  splines  generally  match 
well  at  subdomain  boundaries  but  do  not  always  do  so. 
The  current  hypothesis  is  that  the  cases  in  which  the  local 
Lx  smoothing  splines  do  not  match  well  at  the  boundaries 
of  the  subdomains  are  due  to  limitations  in  the 
compressed  primal-dual  algorithm  that  is  used  to  calculate 
the  local  Lx  smoothing  splines.  The  non-iterative  nature  of 
this  new  domain  decomposition  procedure  is  in  strong 
contrast  to  and  is  a  large  improvement  over  the  iterative 
nature  of  all  previously  known  domain  decomposition 
procedures.  With  sequential  and  especially  with  parallel 
computation,  the  non-iterative  Lx  smoothing  spline 
domain  decomposition  procedure  will  be  a  large  factor  in 
reducing  computing  time  so  that  complex  terrain  models 
can  be  calculated  and  manipulated  in  real  time. 


splines  have  been  used  to  represent  natural  and  urban 
terrain,  geophysical  features  and  financial  processes 
(Champion  and  Lavery,  2002;  Gilsinn  and  Lavery,  2002; 
Lavery,  2001,  2004;  Lavery  and  Gilsinn,  2000,  2001; 
Wang  et  al.,  2006)  and  can  be  used  to  represent  other 
types  of  irregular  data/functions  such  as  geographical 
information,  biological  objects,  mechanical  objects, 
images,  economic  processes  and  social  processes.  Lx 
interpolating  splines  have  dominated  the  previous 
research  on  Lx  splines.  However,  Lx  interpolating  splines 
have  no  compression  capability.  In  contrast,  Lx  smoothing 
splines  do  provide  compression,  at  the  expense  of  some 
loss  of  information,  of  course.  Previously,  Lx  smoothing 
splines  have  been  calculated  globally  (Champion  and 
Lavery,  2002;  Gilsinn  and  Lavery,  2002).  In  this  present 
paper,  we  investigate  a  new,  non-iterative  domain 
decomposition  procedure  for  calculating  Lx  smoothing 
splines  and  present  computational  results  for  urban  and 
natural  terrain.  A  related  non-iterative  domain 
decomposition  procedure  for  Lx  interpolating  splines  is 
currently  under  development  by  the  authors  of  this  present 
paper. 


2.  COMPUTATIONAL  MODEL 

Bivariate  cubic  Lx  smoothing  splines,  proposed  by 
Gilsinn  and  Lavery  (2002),  are  defined  by  minimizing  the 
functional 


1.  INTRODUCTION 

Over  the  past  five  years,  a  number  of  publications 
have  provided  evidence  that  Lx  smoothing  splines  and 
related  Lx  interpolating  splines  preserve  the  shape  of 
urban  and  natural  terrain  well,  better  than  any  other  type 
of  spline  (polynomial,  rational,  exponential  or 
trigonometric).  Cubic  Lx  smoothing  and  interpolating 
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Here,  D  is  a  2D  “global  domain”  over  which  a  piecewise 
cubic  Lx  smoothing  spline  z(x,  y)  is  to  be  calculated.  The 
data  to  be  approximated  are  f  given  at  locations  ( xm ,  j>  ) , 

m  =  1,  2,  M  with  weights  $  .  The  quantity  a  is  a 

balance  parameter  that  determines  the  trade-off  between 
fitting  the  data,  represented  by  the  sum,  and  ensuring  that 
the  smoothing  spline  does  not  have  extraneous  oscillation, 
represented  by  the  integral. 

The  cubic  Lx  smoothing  splines  that  we  use  in  this 
paper  are  created  on  regularly  spaced  rectangular  grids 
with  nodes  (x.,y.),  i  =  0,  1,  ...,  /,  j  =  0,  1,  ...,  J,  on  the 

domain  D  -  (x0,x/)x(y0,yJ)  and  are  based  on  Sibson 

elements  (Han  and  Schumaker,  1997;  Lavery,  2001). 
These  smoothing  splines  are  characterized  by  their 
elevation  z(xi,yj )  and  derivatives  dz(xi,yj)/dx  and 

dz(xi,yj)/dy  at  the  nodes  (xi,yj)  •  The  elevations  and 

derivatives  are  determined  by  minimizing  functional  (1). 
To  carry  out  this  minimization,  we  discretize  functional  (1) 
in  the  following  manner.  Divide  each  cell 
(xi,xi+l)x(yj,yj+l) of  D  into  K2  equal  subrectangles,  K>  2. 

The  integral  over  the  cell  is  approximated  byl/2K(K-l) 
times  the  sum  of  the  2K(K-\)  values  of  the  integrand  at 

the  midpoints  of  the  sides  of  the  subrectangles  that  are  in 
the  interior  of  the  rectangle.  Minimization  of  the 
discretized  version  of  (1)  is  a  linear  program  that  we  solve 
using  a  compressed  primal-dual  algorithm  (Wang  et  al., 
2006).  For  the  computational  results  in  this  present  paper, 
we  added  a  matrix  re-ordering  procedure  to  the 
compressed  primal-dual  algorithm  to  reduce  the 
bandwidth  of  the  reweighted  least  squares  matrices  that 
occur  in  this  algorithm  and  thereby  increase  the  size  of  the 
domains  that  the  algorithm  can  handle. 

Let  I  be  a  divisor  of  /  and  J.  Divide  the  global 
domain  D  into  equal  subdomains  of  LxL  cells.  Except  at 
the  boundary  of  the  global  domain,  extend  each 
subdomain  of  LxL  cells  by  E  cells  in  all  directions  (top, 
bottom,  right  and  left).  This  yields  “extended 
subdomains”  of  (. L  +  2E)x(L  +  2E )  cells  (linear 

dimension  L  +  E  cells  when  a  subdomain  borders  on  the 
boundary  of  the  global  domain).  The  domain 
decomposition  procedure  consists  of  three  steps,  namely, 
1)  calculating  the  Lx  smoothing  spline  individually  on 
each  of  the  extended  subdomains,  2)  restricting  the  spline 
on  each  extended  subdomain  to  the  basic  subdomain  of 
LxL  cells  inside  that  extended  subdomain  and  3) 
creating  the  Lx  smoothing  spline  on  the  global  domain  D 
by  patching  together  the  Lx  smoothing  splines  on  the 
LxL  -cell  subdomains.  In  contrast  to  all  previous  domain 
decomposition  procedures,  this  procedure  is  non-iterative , 
that  is,  it  does  not  require  that  information  be  transferred 
between  subdomains  and  that  Step  1  be  repeated. 
Previous  domain  decomposition  procedures  cannot  be 


non-iterative  because  the  minimization  principles  were 
such  that  a  change  in  the  data  at  any  point  in  the  global 
domain  affects  the  smoothing  spline  everywhere  in  the 
global  domain.  In  contrast,  a  change  in  the  data  at  a  given 
point  affects  the  Lx  smoothing  spline  only  in  a  limited 
region  around  that  point.  Thus,  Lx  smoothing  splines  can 
be  calculated  by  a  non-iterative  domain  decomposition 
procedure  as  long  as  L  and  E  are  sufficiently  large. 
Determining  appropriate  L  and  E  is  a  critical  issue  in 
implementing  an  Lx  smoothing  spline  domain 
decomposition  procedure.  We  consider  this  issue  in  the 
next  section. 

We  will  present  computational  results  for  three 
artificial  data  sets  and  two  real-terrain  data  sets.  The  three 
artificial  data  sets  all  consist  of  641x641  points  on  an 
equally  spaced  rectangular  g  rid.  The  first  artificial  data 
set  has  data  with  height  0  to  the  left  of  a  north-south  line 
through  the  center  of  the  xy  grid  of  this  data  set  and  height 
1  on  and  to  the  right  of  this  line.  The  second  artificial  data 
set  has  data  with  height  0  to  the  left  of  a  northeast- 
southwest  diagonal  line  through  the  center  of  the  xy  grid 
of  this  data  set  and  height  1  on  and  to  the  right  of  this  line. 
The  third  artificial  data  set  has  data  with  height  0  and  1 
alternatively  for  every  10  rows  of  the  xy  grid.  The  first 
real-terrain  data  set  is  an  urban  terrain  data  set  consisting 
of  1000x1000  points  of  lm-spacing  data  for  Baltimore, 
MD  provided  by  the  Joint  Precision  Strike  Demonstration 
Project  Office  (JPSD  PO)  Rapid  Terrain  Visualization 
(RTV)  ACTD.  The  second  real-terrain  data  set  is  a 
natural-terrain  data  set  consisting  of  1201x1201  points  of 
DTED1  (100m  spacing)  data  for  Killeen,  Texas  provided 
by  the  National  Geospatial-Intelligence  Agency.  The 
urban-terrain  and  natural-terrain  data  sets  were  previously 
used  in  (Champion  and  Lavery,  2002;  Lavery,  2001,  2004; 
Lavery  and  Gilsinn,  2000,  2001;  Wang  et  al.,  2006). 

In  all  of  the  computational  experiments  in  this  paper, 
there  are  9x9  =  81  points  in  each  closed  cell 
[x;,x.+]]x[y7.,y7+1]  •  For  very  large  domains,  the  “raw 

compression  ratio”,  that  is  the  number  of  degrees  of 
freedom  in  the  data  divided  by  the  number  of  degrees  of 
freedom  in  the  Lx  smoothing  spline  is  therefore  82/3  = 
21.33.)  For  all  the  computation  results  in  this  paper,  K  =  3 
and  wm- 1  for  m  =  1 ,  2,  . . . ,  M. 


3.  DETERMINING  SIZES  OF  SUBDOMAINS  AND 
EXTENDED  SUBDOMAINS 

3.1  Relation  between  E  and  L 

To  establish  a  relation  between  E  and  L ,  we  will 
calculate  the  E  that  minimizes  the  computational  cost  of 
the  domain  decomposition  procedure  on  a  sequential 
computer.  In  the  compressed  primal-dual  algorithm  that 
we  use  to  calculate  Lx  splines,  the  computational  cost  is 


2 


dominated  by  the  cost  of  the  factorization  of  the 
symmetric  reweighted  least-squares  matrix  on  each 
iteration  of  this  algorithm.  (The  number  of  iterations  of 
the  compressed  primal-dual  algorithm  for  Lx  smoothing 
splines  is  roughly  independent  of  the  size  of  the 
domain/subdomain  and  is  not  a  factor  here.)  Without  loss 
of  generality,  assume  that  J  <1.  For  a  global  domain  of 
IxJ  cells,  there  are  3(7  +  l)x(J  +  l) unknowns  (z(xi,yj) , 

dz(Xi,yj)/dx  and  dz{xt,y  I  dy& each  node  (x i,y])\ the 

number  of  superdiagonals  (and  subdiagonals)  in  the  least- 
squares  matrix  is  3.7  +  9  and  the  factorization  of  this 
reweighted,  least-squares  matrix  therefore  costs,  to 
highest  order,  243T/3  flops.  For  each  extended 
subdomain  of  (. L  +  2E)x(L  +  2E )  cells,  there  are 
3(L  +  2E  +  l)2  unknowns,  the  number  of  superdiagonals  in 
the  least-squares  matrix  of  the  compressed  primal-dual 
algorithm  is  3L  +  6E  +  9  and  the  factorization  of  this 
matrix  costs,  to  highest  order,  243(7, +  2£)4  flops.  The 
total  operation  count  on  a  sequential  computer  for  the 
(7/7,)x(J/7,)  subdomains  is  therefore 

243(7,  +  2E)4  (7  /  L)(J  /  L)  .  With  7,  J  and  L  fixed,  the 
minimum  of  this  expression  occurs  when  E-LI 2.  We 
will  use  this  relation  as  a  basis  for  the  computational 
experiments  reported  in  Subsec.  3.2  below. 

3.2  Determining  E  from  Computational  Experiments 

For  the  domain  decomposition  procedure  outlined  in 
Sec.  2  to  perform  non-iteratively  as  claimed,  the  “width” 
E  of  the  ring  around  each  subdomain  needs  to  be  large 
enough  that  the  Lx  smoothing  spline  on  the  basic 
subdomain  of  LxL  cells  inside  the  extended  subdomain 
is  independent  of  the  data  outside  of  the  extended 
subdomain.  To  determine  how  large  E  needs  to  be,  we 
conduct  a  series  of  computational  experiments  to 
determine  how  far  a  perturbation  in  the  data  propagates  in 
the  Lx  smoothing  spline. 

In  all  of  the  computational  experiments  in  this 
subsection,  the  total  number  of  data  points  is  64 12,  the 
size  of  the  domain  D  is  80x80  cells  (which  corresponds  to 
L  =  40,  E  =  L/2  =  20)  and  there  are  (as  previously  stated) 
9x9  =  81  points  in  each  closed  cell.  For  one  set  of 
computational  experiments,  we  created  a  perturbation  by 
adding  a  constant  p  to  the  heights  at  the  9x9  points  at  the 
center  of  the  data  set.  In  another  set  of  computational 
experiments,  we  added  a  constant  p  to  the  heights  at  a  set 
of  9x9  points  at  a  comer.  To  determine  the  distance  that  a 
perturbation  in  the  data  propagates  in  an  Lx  smoothing 
spline,  we  compare  the  Lx  smoothing  spline  of  the 
perturbed  data  with  the  Lx  smoothing  spline  of  the 
original  data.  Differences  in  heights  at  the  nodes  are 
considered  significant  (“nonzero”)  if  they  are,  in  absolute 
value,  >  10 ~3p.  (The  level  10 ~3p  is  consistent  with  the 


numerical  accuracy  of  the  compressed  primal-dual 
algorithm  used  to  calculate  Lx  smoothing  splines.) 

We  carried  out  3  sets  of  computational  experiments 
with,  respectively,  the  first,  second  and  third  artificial  data 
sets,  each  with  balance  parameter  a  =  0.05,  0.1  and  0.2. 
We  also  carried  out  2  sets  of  computational  experiments 
with  a  641x641  portion  of  the  urban  terrain  data  set  and  2 
sets  of  computational  experiments  with  a  641x641  portion 
of  the  natural  terrain  data  set ,  each  with  a  =  0.05,  0.1  and 
0.2.  In  these  computational  experiments,  the  minimum, 
median,  average  and  maximum  distances  propagated  by 
the  perturbation  in  the  Lx  smoothing  splines  were  0,  4, 
12.19  and  40  cells,  respectively.  When  the  optimal  a  was 
chosen,  the  minimum,  median,  average  and  maximum 
distances  propagated  by  the  perturbation  in  the  Lx 
smoothing  splines  were  0,  2,  5.07  and  11,  respectively. 
(The  optimal  a  value  can  be  different  for  different  data 
sets.)  In  the  following,  we  present  results  of 
computational  experiments  for  the  Baltimore  data  set  with 
p  =  100  and  a  =  0.2  and  for  the  Killeen  data  set  with  p  = 
200  and  a  =  0.2  . 

Figures  la  and  lb  are  the  original  surface  and  the  Lx 
smoothing  spline,  respectively,  for  the  Baltimore  data. 
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Fig.  la.  641x641  portion  of  Baltimore  data  set 
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Fig.  lb.  Lx  smoothing  spline  for  641x641  portion  of 
Baltimore  data  set 
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Figures  2a  and  2b  are  the  original  surface  and  the  Lx 
smoothing  spline,  respectively,  for  the  Baltimore  data 
with  the  perturbation  at  the  center  of  the  data.  Fig.  3  is  the 
difference  between  the  Lx  smoothing  spline  with  the 
perturbation,  shown  in  Fig.  2b,  and  the  Lx  smoothing 
spline  without  the  perturbation,  shown  in  Fig.  lb. 
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Fig.  2a.  641x641  portion  of  Baltimore  data  set  with 
perturbation  at  center 
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Fig.  2b.  Lx  smoothing  spline  for  641x641  portion  of 
Baltimore  data  set  with  perturbation  at  center 


Fig.  3.  Difference  between  the  Lx  smoothing  splines  of 
Figs  2b  and  lb 


Figures  4a  and  4b  are  the  original  surface  and  the  Lx 
smoothing  spline,  respectively,  for  the  Baltimore  data 
with  the  perturbation  at  a  corner  of  the  data.  Fig.  5  is  the 
difference  between  the  Lx  smoothing  spline  with  the 
perturbation,  shown  in  Fig.  4b,  and  the  Lx  smoothing 
spline  without  the  perturbation,  shown  in  Fig.  lb. 
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Fig.  4a.  641x641  portion  of  Baltimore  data  set  with 
perturbation  at  comer 
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Fig.  4b.  Lx  smoothing  spline  for  641x641  portion  of 
Baltimore  data  set  with  perturbation  at  comer 


Fig.  5.  Difference  between  the  Lx  smoothing  splines  of 
Figs  4b  and  lb 
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Figures  6a  and  6b  are  the  original  surface  and  the  Lx 
smoothing  spline,  respectively,  for  the  Killeen  data. 
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Fig.  6a.  641x641 portion  of  Killeen  data  set 
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Fig.  6b.  Lx  smoothing  spline  for  641x641  portion  of 
Killeen  data  set 

Figures  7a  and  7b  are  the  original  surface  and  the  Lx 
smoothing  spline,  respectively,  for  the  Killeen  data  with 
the  perturbation  at  the  center  of  the  data.  Fig.  8  is  the 
difference  between  the  Lx  smoothing  spline  with  the 
perturbation,  shown  in  Fig.  7b  and  the  Lx  smoothing 
spline  without  the  perturbation,  shown  in  Fig.  6b.  Figures 
9a  and  9b  are  the  original  surface  and  the  Lx  smoothing 
spline,  respectively,  for  the  Killeen  data  with  the 
perturbation  at  a  comer  of  the  data.  Fig.  10  is  the 
difference  between  the  Lx  smoothing  spline  with  the 
perturbation,  shown  in  Fig.  9b  and  the  Lx  smoothing 
spline  without  the  perturbation,  shown  in  Fig.  6b. 
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Fig.  7a.  641x641 portion  of  Killeen  data  set  with 
perturbation  at  center 
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Fig.  7b.  Lx  smoothing  spline  for  641x641  portion  of 
Killeen  data  set  with  perturbation  at  center 


3  >.•••■" 


x 


Fig.  8.  Difference  between  the  Lx  smoothing  splines  of 
Figs  7b  and  6b 
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Fig.  9a.  641x641 portion  of  Killeen  data  set  with 
perturbation  at  comer 
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Fig.  9b.  Lx  smoothing  spline  for  641x641  portion  of 
Killeen  data  set  with  perturbation  at  corner 


Fig.  10.  Difference  between  the  Lx  smoothing  splines  of 
Figs  9b  and  6b 

In  Figs.,  3,  5,  8  and  10,  the  significant  differences, 
that  is,  those  >10-3/?  in  absolute  value,  occur  at  most  6 


cells  distant  from  the  locations  of  the  perturbations  in  Figs. 
2a,  4a,  7a  and  9a. 


4.  COMPUTATIONAL  RESULTS  FOR  URBAN 
AND  NATURAL  TERRAIN 

So  that  the  ring  by  which  each  subdomain  is  extended 
is  large  enough  to  assure  that  the  Lx  smoothing  spline  in 
that  subdomain  is  not  affected  by  the  data  outside  the 
extended  subdomain,  we  need  to  choose  E  to  be  greater 
than  the  maximum  of  6  cells  that  the  perturbation 
propagates  in  Figs.  3,  5,  8  and  10  and  greater  than  the 
maximum  of  15  cells  that  the  perturbation  propagates  in 
the  other  computational  experiments  reported  at  the 
beginning  of  Subsec.  3.2.  We  choose  E  =  20.  We  use  the 
optimal  a  =  0.2  . 

We  first  calculate  the  global  Lx  smoothing  spline  for 
a  993x993  portion  of  the  Baltimore  data  set,  depicted  in 
Fig.  11a,  by  domain  decomposition.  We  set  up  a 
smoothing  spline  grid  of  124x124  equal  cells  that 
precisely  covers  the  993x993  data  grid.  We  divide  the 
smoothing  spline  grid  into  four  subdomains,  each  of 
62x62  cells  (Z  =  62),  that  overlap  only  at  their  boundaries. 
We  calculate  the  local  Lx  smoothing  splines  on  each  of  the 
four  extended  subdomains  of  size  82x82  and  then  create 
the  global  Lx  smoothing  spline  by  patching  together  the 
four  local  Lx  smoothing  splines.  Numerical  errors  in  the 
compressed  primal-dual  algorithm  result  in  the  four  local 
Lx  smoothing  splines  not  matching  exactly  on  overlapping 
boundaries.  The  minimum,  median,  average  and 
maximum  of  the  absolute  values  of  the  differences  of  the 
z,  dz/dx  and  dz/dy  at  the  smoothing  spline  nodes  on 
overlapping  subdomain  boundaries  are  given  in  Table  1. 


Table  1.  Minimum,  median,  average  and  maximum  of  the 
absolute  values  of  the  differences  of  quantities  on 
overlapping  subdomain  boundaries  for  Baltimore  data. 


Min 

Median 

Average 

Max 

0 

0.0071 

0.1575 

2.6137 

3  z(xj,yj)/dx 

0 

0.0027 

0.0602 

1.3892 

dz{xl,yj)!dy 

0 

0.0008 

0.0250 

1.5404 

As  the  data  in  Table  1  indicates,  the  differences  in  z, 
dz/dx  and  dz/dy  at  the  smoothing  spline  nodes  on 
overlapping  subdomain  boundaries  are  generally  small 
compared  to  the  total  difference  in  elevation  of  1 00  in  the 
Baltimore  data  set.  The  maximum  absolute  differences  of 
2.6137,  1.3892  and  1.5404  in  z  and  its  derivatives  need  to 
be  investigated  further  to  determine  whether  they  are  due 
to  numerical  limitations  of  the  compressed  primal-dual 
algorithm  or  to  some  other  cause.  When  patching  the  local 
Lx  smoothing  splines  together  to  create  the  global  Lx 
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smoothing  splines,  we  resolve  the  (generally  only  slight) 
differences  in  z,  dz/dx,  and  dz/dy  at  the  nodes  on 
overlapping  subdomain  boundaries  by  averaging  the  two 
or  more  values  of  each  quantity  at  each  node.  In  Fig.  1  lb, 
we  present  the  global  Lx  smoothing  spline  calculated  by 
domain  decomposition  for  the  993x993  portion  of  the 
Baltimore  data  set. 


Fig.  11a.  993x993  portion  of  Baltimore  data  set 


at  the  smoothing  spline  nodes  on  overlapping  subdomain 
boundaries  are  given  in  Table  2. 


Table  2.  Minimum,  median,  average  and  maximum  of  the 
absolute  values  of  the  differences  of  quantities  on 
overlapping  subdomain  boundaries  for  Killeen  data. 


Min 

Median 

Mean 

Max 

0 

0.1600 

0.3556 

4.2500 

dz(xl,yJ )  /  dx 

0 

0.0320 

0.0899 

1.1394 

dz(xi,yj)!dy 

0 

0.0280 

0.0703 

1.4139 

As  was  the  case  for  the  data  in  Table  1,  the  data  in 
Table  1  indicates  that  the  differences  in  z,  dz/dx  and  dz/dy 
at  the  smoothing  spline  nodes  on  overlapping  subdomain 
boundaries  are  generally  small  compared  to  the  total 
difference  in  elevation  of  230  in  the  Killeen  data  set  but 
the  maximum  absolute  differences  of  4.2500,  1.1394  and 
1.4139  need  to  be  investigated  further.  Averaging  the 
(generally  only  slightly)  different  values  of  z,  dz/dx ,  and 
dz/dy  at  the  nodes  on  overlapping  subdomain  boundaries, 
we  produce  the  global  Lx  smoothing  spline  calculated  by 
domain  decomposition  for  the  Killeen  data  set  presented 
in  Fig.  12b. 


decomposition  for  993x993  portion  of  Baltimore  data  set 

Fig.  12a.  1201x1201  Killeen  data  set 


We  also  calculate  the  global  L\  smoothing  spline  for 
the  Killeen  data  set,  depicted  in  Fig.  12a,  by  domain 
decomposition.  We  set  up  a  smoothing  spline  grid  of 
150x150  equal  cells  that  precisely  covers  the  1201x1201 
data  grid.  We  divide  the  smoothing  spline  grid  into  nine 
subdomains,  each  of  50x50  cells  ( L  =  50),  that  overlap 
only  at  their  boundaries.  We  calculate  the  local  Lx 
smoothing  splines  on  each  of  the  four  extended 
subdomains  of  size  70x70  cells  (at  comers),  the  four 
extended  subdomains  of  size  70x90  cells  (on  boundaries 
of  the  global  domain  but  not  at  a  comer)  and  the  one 
extended  subdomain  of  size  90x90  cells  (in  the  interior) 
and  then  create  the  global  Lx  smoothing  spline  by 
patching  together  the  nine  local  Lx  smoothing  splines  on 
the  basic  subdomains  inside  these  extended  subdomains. 
The  minimum,  median,  average  and  maximum  of  the 
absolute  values  of  the  differences  of  the  z,  dz/dx  and  dz/dy 
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Fig.  12b.  Lx  smoothing  spline  calculated  by  domain 
decomposition  for  Baltimore  data  set 


5.  CONCLUSION 

The  computational  results  presented  in  this  paper 
indicate  considerable  success  in  designing  a  new  non¬ 
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